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Numerical  Simulation  of  Three-Dimensional  Self-Gravitating  Flow 


John  V.  Shebalin* 
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The  three-dimensional  flow  of  a  self-gravitating  fluid  is  numerically  simulated  using  a  Fourier  pseu- 
dospectral  method  with  a  logarithmic  variable  formulation  Two  cases  with  zero  total  angular  momentum 
are  studied  in  detail,  a  32^  simulation  (Run  A)  and  a  64®  simulation  (Run  B).  Other  than  the  grid  size,  the 
primary  differences  between  th  ^  two  cases  ate  that  Run  A  modelled  atomic  hydrogen  and  had  considerably 
more  compressible  motion  initially  than  Run  B,  which  modelled  molecular  hydrogen.  The  numerical  results 
indicate  that  gravitational  collapse  can  proceed  in  a  variety  of  ways.  In  the  Run  A,  collapse  led  to  an 
elongated  tube-like  structure,  while  in  the  Run  B,  collapse  led  to  a  flatter,  disk-like  structure. 


'Research  supported  by  the  .National  Aeronautics  and  Spare  Administration.  This  work  was  done  while  the  author  was  in 
residence  as  a  Visiting  Scientist  at  the  Institute  for  Computer  Applications  in  Science  and  Engineering  (ICASE),  NASA  Langley 
Research  Cenicr,  Hampton,  VA  23681. 


1  Introduction 


In  recent  years,  Fourier  methods  have  been  used  to  simulate  collapse  of  a  gas  due  to  a  newtonian  gravitational 
self- interact  ion,  both  in  one-dimension  (ID)  and  in  two-dinieusions  (‘2D).  In  the  ID  case  [1],  the  flow  was 
assumed  to  be  isentropic,  64  grid  points  were  used,  and  it  wa^  found  necessary  to  modify  the  values  of  density 
at  each  time-step  by  adding  a  spatial  constant  so  as  to  ensure  sufficient  positive-definiteness  of  the  density 
at  all  grid  points.  In  the  2D  case  ['i],  the  flow  was  again  assumed  isentropic,  the  maximum  grid  size  was 
‘ioG-.  the  density  was  again  ‘modified’  to  ensure  numerical  stability,  and,  in  addition,  the  flow  was  initially 
turbulent  (with  the  notable  result  of  hinderng  collapse). 

Here,  this  work  is  extended  in  four  ways.  First  the  full  polytropic  gao-dy nai,iic  equrtiioiis  are  solved 
with  no  i.sentropic  approximation.  (Collapse  is  assumed  to  take  place  on  length  and  time  scales  such  that 
no  overall  cosmological  expansion  need  be  imposed  on  the  dynamics.)  Second,  simulations  are  performed  on 
three-dimensional  grids  (32^  and  64'^).  Third,  the  governing  equations  are  written  so  that  the  logarithms  of 
the  density  and  temperature  are  solved  for,  rather  than  the  density  and  temperature  directly.  This  allows 
explicitly  for  the  preservation  of  positive-definiteness  of  density,  temperature,  and  pressure,  since  these  are 
arrived  at  by  exponentiation.  Fourth,  bulk  viscosity  is  utilized  to  model  more  accurately  the  properties  of 
molecular  hydrogen. 

These  extensions  allow  for  more  realistic  simulations  of  the  self-gravitational  collapse  of  clouds  of  either 
atomic  (//)  or  molecular  (Ho)  hydrogen,  at  least  up  to  the  point  where  density  gradients  are  so  steep  that 
the  computer  code  loses  accuracy.  In  particular,  two  cases  with  zero  total  angular  momentum  are  studied 
in  detail,  one  on  a  32^  grid  (Run  A;  H)  and  the  other  on  a  64®  grid  (Run  B;  H2)-  The  primary  difference 
between  the  two  cases  (other  than  grid  size  and  chemical  constituents)  lies  in  the  ratio  of  compressible  kinetic 
energy  to  total  kinetic  energy.  Run  A  had  an  initial  ratio  of  0.5,  while  Run  B  had  an  initial  ratio  of  0.1. 
Collapse  from  a  constant  density  initial  state  is  clearly  different  for  the  two  cases:  in  Run  A,  the  collapse  is 
to  a  tube-like  structure,  while  in  Run  B.  it  is  to  a  flat,  disk-like  structure  (a  Zeldovich  'pancake'  [3])- 

Following  brief  descriptions  of  the  governing  equations  and  numerical  method,  numerical  results  will 
be  discussed.  Flow  parameters  will  be  defined,  as  well  as  measures  of  anisotropy  which  are  useful  for 
quantitatively  describing  self-gravitational  collapse.  Finally,  a  summary  and  conclusion  will  be  presented. 
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Physics  Numerics 


Placing  /<  =  (  '  and  7'  =  f inlo  tin'  standard  polytro|)ic  gits  »’(.)nations  yn-lds  tin-  basic  non-dimensional 
etiuations  in  a  logarithmic  formulation  [1]; 
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Here,  Tij  —  d,Uj  +dju,  — 2/3^,jV  u  and  I  =  is  the  unit  dyadic.  In  (2),  is  the  newtonian  gravitational 
potential  which  satisfies  V'd>  =  k'^p,  where  kj  \s  the  Jeans  wavenumber  [5];  for  Run  A.  fcj  =  \/2  and  for 
Run  B.  kj  —  l.b.  The  adiabatic  index  7  is  constant,  and  since  the  gas  is  either  H  or  7  =  5/3  for  Run  A 
and  7  =  7/5  for  Run  B.  The  dimensionless  shear  viscosity  p  and  thermal  conductivity  k  are  al.so  constant; 
the  Prandtl  number.  Pr  =  7/</k,  has  a  value  of  ~  I,  which  will  determine  k:  k  —  '/p  (p  =  0.02  for  Run  .A 
and  p  =  0.005  for  Run  B).  The  ratio  of  bulk  to  shear  viscosity,  0  =  pb/ p,  "ill  also  be  constant:  =  1  for 
Run  .A  (to  represent  a  small  admixture  of  H2  with  primarily  H)  and  /3  =  32  for  Run  B  (only  H^)  [6]. 

Details  of  the  numerical  method  used  here  have  been  given  previously  [4,  7].  The  only  modification 
lies  in  the  addition  of  the  self-gravitation  term  —  in  (2).  Using  V^<I>  ==  k^p  and  Fourier  transformation 
gives  -Vd>  —  ik'j\ik~'^ p(V.),  which  is  very  straightforward  to  implement  when  (l)-(3)  are  written  in  terms 
of  Fourier  coefficients  (as  is  done  in  the  method  utilized  here).  We  now  proceed  on  to  the  results  of  the 
simulations. 


3  Results 

Hun  A  list'd  a  32*  spatial  grid  and  ran  a  total  of  1530()  timi'-steps  (At'.s),  with  2  cpu-sec/At  on  a  Cray  YMP. 
Run  B  used  a  til*  s[)atial  griil  and  ran  a  total  of  4!)8()()  A/'s,  with  13  cpu-sec  /At,  also  on  a  Cray  YMP. 


Run  A  was  terminated  at  simulation  time  i  =  6.57  when  the  dissipation  wave  number  kp  (as  defined  in  [4]) 
had  increased  so  as  to  equal  the  maximum  wave  number  kmai  =  16  Similarly,  Run  B  was  terminated  at 
simulation  time  t  =  5.U1  when  the  dissipation  wave  number  kp  had  increased  so  as  to  equal  its  maximum 
wave  number  kmar  =  32.  The  microscale  Reynolds  number  R\  (also  defined  in  [4])  ranged  between  17  and 
32  for  Run  A,  and  between  90  and  170  for  Run  B. 

In  the  following,  {Q)  will  denote  the  volume  average  (or  average  value  per  grid  point)  of  any  quantity  Q. 
Thus,  the  various  energies  associated  with  the  collapsing  clouds  of  gas  are  the  kinetic  energy  Ek  =  (pu'/  2), 
the  internal  energy  Ei  =  (pT/[',(7  -  1)]),  the  gravitational  energy  Ea  =  (p<J/2),  and  the  total  energy 
E  =  Ek  +  El  +  Eg-  The  tinre  evolution  of  these  energies  is  shown  in  Figure  1  for  Runs  A  and  B. 

As  is  evident  in  Figure  1,  collapse  did  not  begin  immediately,  but  rather  after  the  initial  turbuleiit  motion 
had  decayed  to  a  lower  level,  in  accord  with  previous  results  [2].  (The  initial  turbulent  velocity  spectra  was 
lu(k)p  ~  k^exp—2k'^/kl,  with  kg  =  i  for  Run  A  and  kg  =  4  for  Run  A.)  However,  once  collapse  had 
progressed  sufficiently,  both  the  mean-squared  dilatation  ^  =  5^(V  and  the  mean-squared  vorticity 

(the  enstrophy)  =  j  ^(V  x  began  to  increase.  This  is  shown  in  Figure  2,  where  it  is  clear  that 
Run  A  had  more  compressible  velocity  initially  than  Run  B  did.  [The  compressible  part  of  the  velocity  is 
Uc(k)  =  t"^kk  •  u(k).]  Note  that  both and  Q  begin  to  increase  near  the  end  of  the  runs. 

The  relative  amount  of  compressible  velocity  in  a  flow'  is  quantified  by  the  ratio  x  =  (“c)  /  (^^)>  whose 
importance  was  first  recognized  by  Passot  and  Pouquet  [8].  The  time  evolution  of  x  for  Runs  A  and  B 
is  shown  in  Fiqure  3.  For  Run  A,  x  —  0.25  before  collapse  initiates,  while  for  Run  B,  x  —  0.02  during 
the  equivalent  time  period.  Although  there  are  many  differences  between  Runs  A  and  B,  the  qualitatively 
different  behavior  of  x  prior  to  collapse  may  be  related  to  the  qualitatively  different  dynamical  behavior  of 
the  gas  in  Runs  A  and  B  during  the  collapse  itself. 

The  results  of  gravitational  collapse  in  Runs  A  and  B  are  seen  in  Figures  4  and  5.  These  figures  represent 
the  densities  of  the  respective  clouds  of  gas  at  the  end  of  each  run  and  are  drawn  by  finding  the  point  of 
maximum  density  in  either  run,  extracting  values  of  density  on  the  three  orthogonal  planes  {x  —  y,  y  —  z.  and 
z  —  x)  which  contain  this  point,  and  then  shifting  the  points  on  these  periodic  planes  so  that  the  maximum 
density  lies  in  their  centers.  In  Figure  4  it  is  clear  that  the  gas  has  collapsed  to  a  tube-like  structure;  in 
Figure  5  it  is  clear  that  the  gas  has  collapsed  to  a  much  flatter,  somewhat  irregular,  disk-like  structure. 
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Ill  order  to  gt't  a  better  picture  of  how  these  structures  developed,  it  is  useful  to  define  the  following 
measures  of  anisotropy  Mj  {j  =  y,  [!.*]: 


I'he  time  evolution  of  the  Mj  for  Runs  A  and  B  are  shown  in  Figure  6.  Note  that  the  differences  in  the  final 
states  of  Runs  A  and  B,  as  shown  in  Figures  4  and  5,  grew  directly  from  their  initial  states,  rather  than 
going  through  a  qualitatively  different  intermediate  stage. 

4  Summary  &  Conclusion 

In  this  work  a  three-dimensional  pseudospectral  computer  code  was  used  to  simulate  the  initial  stages  of 
the  gravitational  collapse  of  clouds  of  either  atomic  or  molecular  hydrogen.  The  full  polytropic  gasdynarnic 
equations  were  solved,  with  the  primary  dynamical  variables  being  Inp,  InT,  and  u.  In  addition  to  shear 
viscosity,  bulk  viscosity  was  also  utilized  to  more  exactly  model  the  properties  of  molecular  hydrogen. 

Two  runs  were  examined  in  detail,  Run  A  (on  a  32^  grid)  and  Run  B  (on  a  64^  grid).  As  the  figures 
show,  the  nature  of  the  collapse  was  different  for  each  run:  in  Run  A,  the  collapse  proceeded  to  a  quasi-lD 
striKture,  while  in  Run  B,  the  collapse  proceeded  to  a  qucisi-2D  structure.  One  primary  difference  between 
the  runs  lay  in  the  relative  amount  of  compressible  turbulent  kinetic  energy  prior  to  collapse.  This  suggests 
that  turbulent  density  fluctuations  may  have  an  important  role  to  play  during  gravitational  collapse. 

There  is  a  large  amount  left  unexplored  in  all  of  this,  with  respect  to  variation  in  initial  conditions, 
transport  properties,  and  grid  size,  or  with  respect  to  the  effects  of  magnetic  fields  or  radiation,  for  example. 
The  sheer  length  of  a  single,  moderately  resolved  simulation  (64®  — ►  ~  200  Cray  YMP  cpu-hours)  precludes 
any  quick  but  thorough  investigation  of  ‘parameter  space.’  Instead,  it  is  clear  that  only  a  slow  and  patient 
effort  (along  with  exponential  increcises  in  computer  power)  will  yield  the  results  and  insight  w'hich  are 
sought. 
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Figure  1.  Energies  vs  time  for  (a)  Run  A  and  (b)  Run  B. 


6 


1  2  ^.3  4  5  6 

Tune 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  No  0704  0188 


Pubiii.  ^eportm^  burden  for  collection  o#  mformdtion  is  estimated  to  a«eraqe  ’  nour  per  fesporse.  including  tee  t  me  tor  rewiewng  mstruaions  searcn.ng  anviinq  oata  sources 
gathering  and  mamtainirsg  the  data  needed,  and  comoieting  and  reviev.  ng  the  :oHectionot  mtormation  Seng  comments  reoarding  this  burden  estimate  or  an*  jther  aspect  o’  this 
collection  ot  information,  including  suggestions  for  reducing  this  ouroen  to  wvash>ngion  Heaoc  carters  Services.  Oireaorate  for  information  Operations  ano  Reports,  li  iS  ;eHefson 
OavisHigi^wav.  Suite  U04.  Arlington.  vA  22202-4302.  and  to  the  Office  o»  Management  arsd  Budget.  Paperwom  Reduction  Project  (Of04-OiM).  Washington  DC  ?0S03 


1.  AGENCY  USE  ONLY  (Leave  blank) 


2.  REPORT  DATE 

3  REPORT  TYPE  ANO  DATES  COVERED 

August  1993 

Technical  Memorandum 

4.  TITLE  ANO  SUBTITLE 

NUMERICAL  SIMULATION  OF  THREE-DIMENSIONAL  SELF-GRAVITATING 
FLOW 


5.  FUNDING  NUMBERS 


6.  AUTHOR(S) 

John  V.  Sheballn 


WU  505-90-52-01 


7.  PERFORMING  ORGANIZATION  NAME(S)  ANO  AOORESS(ES) 
Institute  for  Computer  Applications  in  Science 
and  Engineering 

Mail  Stop  132C,  NASA  Langley  Research  Center 
Hampton,  VA  23681-0001 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

ICASE  Report  No.  93-60 


9.  SPONSORING /MONITORING  AGENCY  NAME(S)  AND  AODRESS(ES) 
National  Aeronautics  and  Space  Administration 
Langley  Research  Center 
Hampton,  VA  23681-0001 


10.  SPONSORING /MONITORING 
AGENCY  REPORT  NUMBER 

NASA  TM-109021 
ICASE  Report  No.  93-60 


11.  SUPPLEMENTARY  NOTES 

Langley  Technical  Monitor:  Michael  F.  Card 
Final  Report 


To  be  submitted  to  Astronomy  and 
Astrophysics 


IZa.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 

Unclassified  -  Unlimited 


12b.  DISTRIBUTION  CODE 


Subject  Category  3A 


13.  ABSTRACT  (Maximum  200  words) 

The  three-dimensional  flow  of  a  self-gravitating  fluid  is  numerically  simulated. using 
a  Fourier  pseudospectral  method  with  a  logarithmic  variable  formulation.  Two  cases 
with  zero  total  angular  momentum  are  studied  in  detail,  a  32^  simulation  (Run  B) . 
Other  than  the  grid  size,  the  primary  difference  between  the  two  cases  are  that  Run  A 
modelled  atomic  hydrogen  and  had  considerably  more  compressible  motion  initially  than 
Run  B,  which  modelled  molecular  hydrogen.  The  numerical  results  indicate  that  gravi¬ 
tational  collapse  can  proceed  in  a  variety  of  ways.  In  the  Run  A,  collapse  led  to  an 
elongated  tube-like  structure,  while  in  the  Run  B,  collapse  led  to  a  flatter,  disk¬ 
like  structure. 


14.  SUBJECT  TERMS 

self-gravitating  flow,  polytrope,  spectral  method 


17.  SECURITY  CLASSIFICATION 
OF  REPORT 

Unclassified 


NSN  7540.01-280-5500 


18.  SECURITY  CLASSIFICATION 
OF  THIS  PAGE 
Unclassified 


19.  security  classification 
OF  ABSTRACT 


IS.  NUMBER  OF  PAGES 

lA 


16.  PRICE  CODE 

A03 


20.  LIMITATION  OF  ABSTRACT 


Standard  Form  298  (Rev  2-89) 

Prmcribed  by  ANSi  Std  239-i6 
298- 102 


☆  l!.S  government  printing  OEFICE:  1993  •  528-064/86070 


